/*It computes the counterfactual density distribution following the parametric method*/
capture program drop density_parametric_compute 
program define  density_parametric_compute 

global data_in=20
global data_out=1180  
global bin=10
global binname=$bin 
global zoom_in_c300=80        
global zoom_out_c300=470
global range_c300   "incbin${binname}>=float(${zoom_in_c300}) & incbin${binname}<=float(${zoom_out_c300})" 



clear mata
matrix table=J(2,8,.)
global rr=1

use  "${density_data_dir}\density_${firm_name}_${name}_${year_name}_bin${binname}_in${data_in}_out${data_out}.dta", clear

global notch_s= ${cutoff}/$bin  
global low_s=${low}/$bin          //normalized;
global high_s=${high}/$bin  

global low_name    =-${low}
local  incbin_name ="incbin${binname}"
local  freq_name   ="numbin${binname}"   //original data

gen round50=0
replace round50=1  if mod(incbin${binname},50)==0

gen Z100=0
replace Z100=1     if mod(incbin${binname},100)==0   
		
keep if $range_c300


/* MATA ESTIMATION */		
mata:mata clear
mata:incbinname=st_local("incbin_name")   
mata:freqname=st_local("freq_name")    
mata:notch_s=strtoreal(st_global("notch_s"))  
mata:low=${low_s}
mata:high=${high_s}
mata:degree=${degree}
mata:nboot=${nboot}
mata:sf_s=$bin  
set seed 1234

* CALLING MATA PROGRAM TO ESTIMATE THE COUNTERFACTUAL, b, a*, z_u AND THE ELASTICITIES *
if "${density}" =="DP_q" & "${firm}"=="domestic"	{ 	
mata: notch_ct_roundq(incbinname,freqname,notch_s,low,high,degree,nboot,sf_s)
}


/* POPULATING THE TABLE */
matrix table[1,1]=.
matrix table[1,2]=${cutoff}
matrix table[1,3]=$bin  
matrix table[1,4]=${degree}
matrix table[1,5]=${low}   //before normalize
matrix table[1,6]=${high}
matrix table[1,7]=${b}
matrix table[2,7]=${b_se}
matrix table[1,8]=${m}
matrix table[2,8]=${m_se}
	
noi di "b=" ${b}
noi di "b_se=" ${b_se}
noi di "m=" ${m}
noi di "m_se=" ${m_se}


/* GETTING THE BUNCHING PLOT */
svmat  y_actual_noroundsround 
svmat  y_counter
svmat  y_counter_noroundsround  

ren y_actual_noroundsround   numbin${binname}_noround          
ren y_counter_noroundsround  numbin${binname}_counter_noround  
ren y_counter                numbin${binname}_counter                 

save "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", replace 
     
svmat table
keep table* 	
rename table1  Region
rename table2  cutoff
rename table3  bin
rename table4  degree 
rename table5  low
rename table6  high
rename table7  b
rename table8  m
	
tostring Region,replace
replace Region = " " if Region== "." 
replace Region="${name}"  in 1  
save "${density_data_dir}/Table1_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${bin}_in${zoom_in_c300}_out${zoom_out_c300}.dta",replace 
		

		
/*draw figs*/
global rr=1	
global b=strofreal(b[${rr}],"%9.2f")	
global b_se=strofreal(b[${rr}+1],"%9.2f")
global m=strofreal(m[${rr}],"%9.2f")	
global m_se=strofreal(m[${rr}+1],"%9.2f")
global graph_low=${low}+${cutoff}
global graph_high=${high}+${cutoff}

end 


capture program drop density_parametric_fig 
program define density_parametric_fig

use "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", clear


#delimit;
twoway 
(connected numbin${binname}                     incbin${binname}, sort clcolor(green) mcolor(green) msymbol(o))
(line      numbin${binname}_counter             incbin${binname}, sort clcolor(maroon)  mcolor(maroon) msymbol(o) lwidth(medthick) lpattern(dash))
,
legend(label(1 "Obesrved")   label(2 "Counterfactual") size(small))
xline(${cutoff}, lcolor(black) lwidth(medthick)) 
xline(${graph_low}  ${graph_high}, lcolor(black) lwidth(medthin)  lpattern(dash)) 
xlabel(${zoom_in_c300} (20)  ${zoom_out_c300}, labsize(small))
ytitle(Density, size(small) )  	
xtitle(Registered Capital(unit:10,000 RMB), size(medsmall) margin(medium) ) xtitle(, alignment(top)) 						
yscale(r(0))	ylabel(#5) 
title("Density distribution of registered capital for", justification(center) color(black) size(medsmall)) 
subtitle( "${firm_name} firms in ${location} ${year_note}"  " ${cutoff_note}",  justification(center) color(black) size(med))
note("bin=${bin}, degree=${degree}, b=${b}(${b_se}).",size(small) )
graphregion(fcolor(white)  ifcolor(white) color(white) icolor(white))
;						
graph export 
"${density_graph_dir}/fitfig_density_${density}_${cutoff}_${firm_name}_${name}_${year_name}.pdf", replace;
#delimit cr		





//no rounding
global zoom_in_c300_rv=160       
global zoom_out_c300_rv=470
global range_c300_rv   "incbin${binname}>=float(${zoom_in_c300_rv}) & incbin${binname}<=float(${zoom_out_c300_rv})"  
#delimit;
twoway 
(connected numbin${binname}_noround                    incbin${binname}, sort clcolor(green) mcolor(green) msymbol(o) msize(large))
(line      numbin${binname}_counter_noround           incbin${binname}, sort clcolor(maroon)  mcolor(maroon) msymbol(o) lwidth(medthick) lpattern(dash))
if $range_c300_rv
,
legend(label(1 "Observed (no rounding)")   label(2 "Counterfactual (no rounding)")  size(small))
xline(${cutoff}, lcolor(black) lwidth(med)) 
xline(${graph_low}  ${graph_high}, lcolor(black) lwidth(medthin)  lpattern(dash)) 
ytitle(Density, size(small)   margin(medium))  
xtitle(Registered Capital(unit:10,000 RMB), size(medsmall))xtitle(, alignment(top)) 						
xlabel(${zoom_in_c300_rv}  (20)  ${zoom_out_c300_rv}, labsize(small))
yscale(r(0))	ylabel(#5) 
graphregion(fcolor(white)  ifcolor(white) color(white) icolor(white))
title("Density distribution of registered capital for", justification(center) color(black) size(medsmall)) 
subtitle( "${firm_name} firms in ${location} ${year_note}"  " ${cutoff_note}",  justification(center) color(black) size(med))
note("Note: bin=${bin}, degree=${degree}, b=${b}(${b_se}).",  size(small))
;
graph export 
"${density_graph_dir}/fitfig_density_${density}_norounding_${cutoff}_${firm_name}_${name}_${year_name}.pdf", replace;
#delimit cr

end 



//exclusion region only
capture program drop density_parametric_exfig
program define density_parametric_exfig

use "${density_data_dir}\density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_degree${degree}_low${low}_high${high}_bin${binname}_in${zoom_in_c300}_out${zoom_out_c300}.dta", clear


#delimit;
twoway 
(connected numbin${binname}                    incbin${binname}, sort clcolor(green) mcolor(green) msymbol(o) msize(large))
(line      numbin${binname}_counter          incbin${binname}, sort clcolor(maroon)  mcolor(maroon) msymbol(o) lwidth(medthick) lpattern(dash))
if incbin${binname}>=float(${graph_low}) & incbin${binname}<=float(${graph_high})
,
legend(label(1 "Observed")   label(2 "Counterfactual")  size(small))
xline(${cutoff}, lcolor(black) lwidth(med)) 
ytitle(Density, size(small)   margin(medium))  
xtitle(Registered Capital(unit:10,000 RMB), size(medsmall))xtitle(, alignment(top)) 						
xlabel(${graph_low} (${binname})  ${graph_high}, labsize(small))
yscale(r(0))	ylabel(#5) 
graphregion(fcolor(white)  ifcolor(white) color(white) icolor(white))
title("Density distribution of registered capital for", justification(center) color(black) size(medsmall)) 
subtitle( "${firm_name} firms in ${location} ${year_note}"  " ${cutoff_note}",  justification(center) color(black) size(med))
note("Note: bin=${bin}, degree=${degree}, b=${b}(${b_se}).",  size(small))
;
graph export 
"${density_graph_dir}/fitfig_density_${density}_${cutoff}_${firm_name}_${name}_${year_name}_excluded.pdf", replace;
#delimit cr
end  